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ABSTRACT 

We perform three dimensional radiation hydrodynamic simulations of the structure and dynamics 
of radiation dominated envelopes of massive stars at the location of the iron opacity peak. One 
dimensional hydrostatic calculations predict an unstable density inversion at this location, whereas our 
simulations reveal a complex interplay of convective and radiative transport whose behavior depends 
on the ratio of the photon diffusion time to the dynamical time. The latter is set by the ratio of 
the optical depth per pressure scale height, tq, to Tc = c/cg, where Cg « 50 km s“^ is the isothermal 
sound speed in the gas alone. When tq ^ Tc, convection reduces the radiation acceleration and 
removes the density inversion. The turbulent energy transport in the simulations agrees with mixing 
length theory and provides its first numerical calibration in the radiation dominated regime. When 
To ^ Tc, convection becomes inefficient and the turbulent energy transport is negligible. The turbulent 
velocities exceed Cg, driving shocks and large density fluctuations that allow photons to preferentially 
diffuse out through low-density regions. However, the effective radiation acceleration is still larger 
than the gravitational acceleration so that the time average density profile contains a modest density 
inversion. In addition, the simulated envelope undergoes large-scale oscillations with periods of a 
few hours. The turbulent velocity field may affect the broadening of spectral lines and therefore 
stellar rotation measurements in massive stars, while the time variable outer atmosphere could lead 
to variations in their mass loss and stellar radius. 

Subject headings: stars: massive — hydrodynamics — methods: numerical — radiative transfer 


1. INTRODUCTION 

Massive stars play an important role in many astro- 
physical environments. The radiative and mechanical 
energy output from massive stars are important feedback 
mechanisms regulating star formation and the structure 


of the inters tellar medium in galaxies (Kennicutt 2005 
Smith 2014). Radiation from massive stars also plays 
an important role in th e reionization of the universe 


(Bromm & Larson 2004). When massive stars explode, 
they produce various types of supernovae and high en¬ 
ergy transients, while leaving behind black holes and neu¬ 
tron stars whose properties depe nd on the previou s his¬ 
tory of massive stellar evolution (jHeger et al.||2003|). 

Because of the steep mass-luminosity relation of 
main sequence stars, the luminosity approaches the 


Eddington-limit for electron scatterir 
masses larger than ~ 50 — lOOMm 

ig in stars 
(Crowther 

with 

2007 

Maeder et al. 

2012 

Sanyal et al. 

2015) 

In massive star 

envelopes witJ 

1 non-zero metallici 

y, the opacity increases 


significantly compared with electron scattering due to 
iron opacit y at temperatures « 1.8 x 10® K (e.g., Paxton 
et al. 2013|and references therein). As a result, the actual 
radiation force can locally exceed the gravitational force 
near the iron opacity peak. This same physics can oper¬ 
ate at somewhat lower temperatures d ue to the opac ity 
produced by helium and hydrogen (e.g. Maeder|fl981 ). 


The locally super-Eddington luminosity near the iron 
opacity peak (and related opacity peaks) has been a long 


standing challenge for one dimensional (ID) stellar evo¬ 
lution models of massive stars. Hydrostatic and ther¬ 
mal equilibrium models with a super-Eddington l umi- 


nosity require density and gas pressure inversions (Joss 


et al. |1973al IGrM'ener et ai.||2012t [Paxton et ai.||2U13t 

Owocki 2015), which are usually convectiv 

ely unstable, 
deep in the 
ge the stellar 
ing the den- 
rium models, 
surface, dow¬ 
ry the neces- 
w the stellar 
ficular if the 

When tne iron opacity peak is sufficiently 
star, convection is efficient and can rearran 
structure to have constant entropy, remov 
sity inversion predicted by radiative equilib 
When the iron opacity peak is close to the 
ever, convection is inefficient and cannot car 
sary super-Eddington flux. It is unclear ho 
structure adjusts in this limit, and in par 

density inversions are stable (Maeder||2009 

Sanyal et al. 

20151. 'Phis results in very uncertain stellar radii in the 

very massive star regime (Yusof et al.|2013 

Kohler et al. 

2015). To complicate the situation even tu 
evolution calculations of massive stars bee 
cally difficult when the opacity near the s 
sphere significantly exceeds the electron sea 
ity. In this situation the convective energy 
often artificially enhanced in ID stellar ev( 

rther, stellar 
ome numeri- 
tellar photo- 
Tering opac- 
transport is 
rlution mod- 
leads to the 

els (Maeder|1987 Paxton et al.|2013), whici 


‘Einstein Fellow 


1979). 

The density and temperature dependence of the opac¬ 
ity near composition-specific opacity peaks can also drive 
radiatio n hydrodynamic instab i lities distinct from con- 
vection (Kiriakidis et al. 1993 Blaes & Socrates 2003 
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Suarez-Madrigal et al. 2013). In many cases, these can- 
not be captured by ID stellar evolution models. These 
instabilities have been suggested as the origin of the large 
mass loss rates from massiv e stars that cannot be ex - 
plained by line-driven winds (Glatzel &: Kiria kidis|1993l). 
The a mount of mass loss by line driven winds (|Fuls et al. 
20081 also depends on the structure of the stellar photo¬ 


sphere (e.g., the amplitude of density and velocity fluc¬ 
tuations, potentially affecting the degree of wind clump- 

ing). 

In this paper, we present three-dimensional (3D) radia¬ 
tion hydrodynamic simulations of stellar envelopes near 
the iron opacity peak and quantify the impact of con¬ 
vective and radiation hydrodynamic instabilities on the 
stellar envelopes’ properties and dynamics. These cal¬ 
culations provide the first di rect calibr ation of classical 
mixing length theory (MLT) (|Cox|1968|) in the radiation 
pressure dominated regime. In addition, we quantify 
the extent to which convection and/or other hydrody¬ 
namic instabilities generate a “porous” atmosphere that 
substantially increases the radiation flux relative to that 
predicted by one-dimensional models. This is critical for 
determining whether locally super-Eddington regions in 
a stellar envelope can initiate a large-scale stellar wind 
via continuum radiation pressure or whether the stel- 
lar radiation freely escapes through low density channels 
( Shaviv|[i998 van Marie et al.|[2008 ). 

Uur work utilizes numerical techniques for accurate 3D 


Eddington tensor (VET) formalism 1 

Jiang et al. 

2012 

Davis et al.||2012|). These techniques 

have been success- 


tulfy used to study a wide range o f astrophysical prob¬ 


lems (Jiang et al. 2013albl 2014bI and have significant 


advantages relative to ad-hoc closures such as the flux- 


limited diffusion (F ED) approximation (Davis et al. 


2014 


Proga et al. 2014). The direct solution of the radiation 


transfer equation (via VET or Monte-Carlo methods) 
is particularly important for studying mildly optically 
thick regions near the stellar photosphere. Indeed, recent 
studies have shown that approximate radiation transfer 
schemes such as FED or Ml can reach completely op¬ 
posite co nclusions compared with VET or Monte Carl o 
methods (Davis et aH2014 Tsang & Milosavljevic|2015 ). 

The remainder of this paper is organized as follows. 
In we identify the key dimensionless parameter that 
governs the physics of near surface convection zones in 
stellar envelopes. We also describe the ID stellar models 
that motivate the specihc conditions in our 3D radiation 
hydrodynamic simulations. In ^we describe our numer¬ 
ical methods while in (|4 we present our results. Finally, 
in (J^we summarize ouiTcey results and discuss their im¬ 
plications. 


2. THE MODEL PARAMETERS 

The relative importance of convective and radiative en¬ 
ergy transport can be characterized by a critical optical 
depth Tc over the pressure scale height Hq near the iron 
opacity peak. As we will show based on our simulations, 
for conditions typical of massive stars, the maximum tur¬ 
bulent velocity reaches the isothermal sound speed when 
convection starts to become inefficient. Then the criti¬ 
cal optical depth can be estimated by balancing the dif¬ 
fusive radiation flux ~ OrTpc/r and the maximal con¬ 
vective flux in the radiation pressure dominated regime 




'7^4. 

. 


Tc — cjCg^Q. 


( 1 ) 


We define the gas isothermal sound speed Cg^, radiation 
sound speed Cs,o, pressure scale height Hq, optical depth 
To and typical dynamic time scale to based on the density 
Po and temperature Tq at the iron opacity peak 



clo arT^ 

^0 — - — n -’ 

9 5po9 

To = Kt{po, To)poHo, 

to = —. (2) 

^s,0 


The characteristic radiation pressure is P^.o = cirT^/S 
while the gas pressure is Pq. 



log (Teff/K) 

Figure 1. Stellar evolution tracks of two solar metallicity stars 
with initial mass of 4OM0 and SOMq calculated with the ID stellar 
evolution code MESA. The stellar symbols identify the locations 
corresponding to the three initial conditions of the 3D radiation 
hydro calculations with Athena listed in Table The Athena 
calculations are effectively local plane-parallel atmosphere calcula¬ 
tions motivated by the MESA models. 


A related way to interpret the critical optical depth Tc 
in equation Q is that the typical thermal time scale due 
to photon dinusion near the iron opacity peak is given 
by 


tc — - —-to- 

C Tq Cp,0 


( 3 ) 


When To ^ Tc, the diffusion time scale is longer than 
the local dynamic time scale and we expect convection 
to be an efficient energy transport mechanism. When 
To <C Tc, however, radiative diffusion is rapid enough 
that convection may behave very differently compared 
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with classical mixing length theory. This can also be 
related to the Peclet number Pe, which is the ratio of 
advective to diffusive flux. When tq ^ Tc, we expect 
Pe > 1, while when tq <C Tc, Pe 1. Note that there 
is some ambiguity as to whether should be defined 
in terms of the gas isothermal sound speed Cg^ or the 
radiation sound speed Cs,o- As we shall show, we find that 
the former better characterizes the maximum velocities 
reached in our simulations. 

In order to identify stellar conditions spanning the dif¬ 
ferent regimes of radiation dominated convection, we use 
Modules for Experiments in St ellar Evolution (MESA, 


Paxton et al. (2011 2013 2015) release r7385) to evolve 
stars with initial niass AOMq and 80Mq. The models 
have an initial metal l icity Z = 0.02 with a mixture taken 


from Asplund et al. ( 2005) and are calculated using the 
OPAL opacity tables (jlglesias & Kogers||iyyt)l). Convec¬ 
tion is accounted for usin g mixing-length theory (MET) 
in the Henyey et al. (19651 formulation with omlt = 1-5. 
Convective boundaries are determined using the Ledoux 
criterion an d we use semiconvection with an efficiency 
Osr = 0.01 (iLanser et al.||1983| 

1|h 


decaying overshooting 


miD 


19851 ). An exponentially 
has been included 
convective bound- 


lerwig 

with ttov = 0.001 above and below 
aries. Stellar winds are i ncluded according to t he mass 
loss scheme described in Glebbeek et al. (2009) (Dutch 
wind scheme) with an efnciency parameter r] = 0.8. 

We have identified specific times in the MESA massive 
star models with conditions that range from tq/tc 1 to 
To/Tf, 3> 1. We label these models as StarDeep, StarMid 
and StarTop; they have tq/tc = 15.2,0.97 and 0.025, re¬ 
spectively. EigureH shows the positions of these models 
in the HR diagramMVIodel StarDeep has an initial mass of 
4OM0. It is evolved to an age of 4.3 x 10® year in MESA 
at which point it has an effective temperature 7.13 x 10® 
K, radius dOlP© and luminosity 8.96 x lO^L©. Model 
StarMid starts from a 80Mq star. At an age of 3.02 x 10® 
year it has an effective temperature 1.17 x 10^ K, ra¬ 
dius 274P© and luminosity 1.26 x 10®L©. Model StarTop 
also starts from an 80M© star but is only evolved to an 
age of 9.15 x 10'^ year, at which point it has an effective 
temperature 4.75 x 10^ K, radius 14P© and luminosity 
8.96 X 10®L©. 

As described in we use plane-parallel 3D radiation 
hydrodynamics simmations to study the physics of mas¬ 
sive stellar envelopes. The initial radial location of the 
iron opacity peak tq and the density po, temperature Tq 
and luminosity Lq at this location are taken from the 
corresponding MESA models. The fiducial parameters 
of our initial conditions are listed in Table [T] 


3. NUMERICAL METHOD 
3.1. Equations 

We model the envelopes of massive stars in plane paral¬ 
lel geometry, neglecting the global stellar geometry. This 
is suitable for studying the effects of 3D radiation hy¬ 
drodynamics on energy transport in the stellar envelope, 
but is not appropriate for studying global effects such 
as stellar outflows. We solve the frequency-independent 
(gray) radiation hydrodynamic equations in cartesian co¬ 
ordinates (x,y,z) with unit vectors (i, y, z). The grav¬ 
itational acceleration g is assumed to be a constant and 
along the — z direction. As the mass of the envelope is less 


than 1% of the core mass, this is a reasonable assumption 
(the slight change in l/r^ in the stell ar envelope is also 


neglec ted). The equations are (e.g., Jiang et al. 2012 
mSal) 


djpv) 

dt 

dE 


f+ v.M = i), 

-I- V • {pvv + PI) = -Sr{P) - pgz, 


— + ^ -[[E + P)v] = -cSr{E) - gpv ■ z, 
BE 

■Fr=cSr{E), 


1 BFr 


■V-P, = Sr(P), 


dt 

where the radiation source terms are 


(4) 


Sr(P) = -p [naF + Ksf) [Fr - {vEr -|- V ■ P^)] /c 
+ pv{KaParT'^ - KaEEr)/c, 

Sr{E) = p{KaparT'^ - KapEr) 

4” Pi^aF ^sF^^ 2 ‘ {vEj. -\- V • P^.)] . (5) 

Notice that local thermal equilibrium (LTE) is assumed 
for the emission term. 

In the above equations, p,P,v,c are the gas density, 
pressure, flow velocity and speed of light respectively. 
The total gas energy density is E = Eg + pv"^ 12, where 
Eg = P /(7 — 1 ) is the internal gas energy density with 
a constant adiabatic index 7 = 5/3. The gas pressure 
is P = pk^T/p, where ke is Boltzmann’s constant and 
p = 0.62mp is the mean molecular weight for nearly fully 
ionized gas with proton mass rUp. The radiation constant 
is Or = 7.57 X 10®® erg cm“® K“^, while Er,Fr are the 
radiation energy density and flux. The Rosseland mean 
absorption and scattering opacities are denoted by KaF 
and KsF, while KaP and KaE are the Planck and energy 
mean absorption opacities. We describe our choice of 
these opacities in the following section. 

Equations (|^ and ([^ are solved in the mixed frame, 
meaning that the radiation flux Fr and energy density 
Er are Eulerian variables, while the material-radiation 


interaction terms in 

equations ([5 

) are written in the co- 

moving frame (e.g., 
and co-moving flux i 

Lowrie et a 


1999 

). The Eulerian 

^^^0 are relatec 

1 

II ■ 

0 

0 


{vEr + V ■ Pj.). Notice that only the co-moving radiation 
flux contributes to the radiation acceleration while the 
advection flux does not. 

We relate the radiation pressure P,. and radiation en¬ 
ergy density through a variable Eddington tensor (VET), 
P^ = fP,,. We do not use an assumed closure relation 
(such as ELD or Ml) to specify f. Instead, we compute 
it directly from angular quadratures of the specific in¬ 
tensity Ir for each angle n, which is calculated from the 
time-independent radiation transfer equation 

BT 

^=^t{S-Ir), (6) 

where S is the radiation source term and Kt is the to¬ 
tal specific opacity. Because the typical dynamic time 
scale is much longer than the light crossing time, f barely 
changes between each time step. Our use of the time 
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independent transfer equation to calculate the angular 
distribution of the radiation field is thus ve^ accurate. 
We solve the transfer equation in equation Q usi ng the 
method of s hort characteristics, as described in |Davis| 
et al.| (|2012|). Then f is calculated according to its deh- 
nitioii 


f = 


f IrUndil 
Jirdn ’ 


( 7 ) 


where is the solid angle. We emphasize that we only 
use the time-independent transfer equation ([6 1 to calcu¬ 
late the Eddington tensor f, which represents the angular 
distribution of the radiation field. The moments of are 
not used to calculate the radiation-material interactions. 
Instead, the radiation moments Er, and the radiation 
source terms are calculated self-consistently through the 
time dependent radiation moment equations Q. 

The Eddington tensor in equation 0 closes the radi¬ 
ation moment equations (eqs. We solve these equa¬ 
tions using the Godunov radi ation MHD c o de At hena, 
as described and tested in |Jiang et al.| (|2012|) (and 
with ad ditional improvements described in | Jiang et al.| 


(2013c)). We solve the radiation subsystems as well as 
the radiation source terms implicitly, while the hydro 
equations are solved explicitly as in the original Athena 


code (Stone et al. 


due to the stitt rac 


2008 ), with approp riate modification s 
lation source terms (Jiang et al.|2012). 


3.2. The Opacity 


In order to correctly capture the temperature and den¬ 
sity dependencies of the opacity, particularly near the 
iron op acity peak, we directly use the op acity table from 
MESA (|Paxton et al.||2011| |2013| |2015 ) in our Athena 
simulations. We assume a rnetalficity Zj — 0.02 and hy¬ 
drogen fraction X = 0.6. For the density p and tem¬ 
perature T in each cell, the total specific opacity kj is 
calculated with bilinear interpolation based on the ta¬ 
ble. This is used as the total specihc flux mean opacity 
Nt = NaF + ^sF- In Order to split this into scattering 
and absorption opacities, we assume a constant electron 
scattering opacity Ksf = 0.32 g“^ cm^. All of the ab¬ 
sorption opacities (koF, FaP, Fas) are assumed to be the 
same for simplicity. 

Figurej^shows the actual opacity k* we use for a range 
of temperatures and densities appropriate to stellar en¬ 
velopes. The opacity peak between 10® and 2 x 10® K 
is due to Fe, which is the region we focus on. The in¬ 
crease in opacity at lower temperatures is due to He and 
H. The iron opacity peak has a sensitive dependence on 
temperature but a weak dependence on density. Partic¬ 
ularly when p < 3.6 X 10“^^ g cm“®, does not de¬ 
pend on density for T > 10® K. For the stellar envelopes 
we study, the temperature usually does not drop below 
4 X 10^ K. Moreover, when this does happen, the den¬ 
sity is always smaller than 10“^® g/cm®. As a result, the 
He and H opacity peaks usually do not show up in our 
calculations. 


3.3. Initial and Boundary Conditions 

We set our initial conditions in the Athena simula¬ 
tions by using the MESA models described in ^to de¬ 
termine the initial radial location rg , density po and tem¬ 
perature To at the iron opacity peak. The gravitational 



Figure 2. The opacity Kt ^ function of density and temperature 
near the iron opacity peak at T ps 1.5 x 10^ K. Each line shows 
the opacity as a function of temperature for a fixed density. From 
top to bottom, the density decreases by a factor of 10 between 
neighboring lines. 


acceleration g is also chosen to be the value at rg in the 
MESA models. We then apply a constant radiation flux 
Fr^i through the whole simulation box, which represents 
the radiation coming from the center of the massive star. 
The initial radiation flux is determined by the total lumi¬ 
nosity Lq in the MESA models at tq: Fr^i = To/(47rrg). 

With these hducial parameters and the opacity given 
in §3.2[ we construct our initial conditions by solving for 
a hydrostatic envelope in thermal equilibrium. Starting 
from z = rg, all the fluid and radiation quantities are 
integrated as a function of height using 


dP 

dz 


PKt p 

=- Fr,i - pg, 

c 


Er = UrT'^, 


1 dEr _ pKt ^ 
3 dz ~ c 


( 8 ) 


We integrate above and below rg to capture the whole 
profile of the iron opacity peak region in our simulation 
box. We then add random perturbations with 0.01% 
amplitude to the density in the simulation box of dimen¬ 
sions Lx, Ly, Lz. The initial conditions constructed in 
this way are slightly different from the ID profiles re¬ 
turned by MESA models around the iron opacity peak. 
We only guarantee that the characteristic physical pa¬ 
rameters we choose are consistent with the MESA stel¬ 
lar models. In particular, our initial condition does not 
have any mixing length model of convective or advective 
energy transport. We then calculate how the energy will 
be transported using self-consistent 3D radiation hydro- 
dynamic simulations. 

At the bottom of the box, we use a reflecting boundary 
condition for gas quantities. This means that the density, 
opacity and the horizontal components of the velocity are 
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Table 1 

Simulation Parameters 


Variables/Units 

StarDeep 

StarMid 

StarTop 

tq/Rq 

125.7 

65.6 

13.6 

To/K 

1.87 X 10® 

1.67 X 10® 

1.57 X 10® 

Po/(g cm-3) 

3.10 X 10“® 

3.60 X 10“® 

5.52 X 10“^ 

g/(cra s'^) 

63.59 

3.56 X 10® 

1.17 X 10® 

FV,i/(erg cm“2 s“^) 

1.24 X 10^2 

9.94 X 10®® 

3.06 X 10®® 

f/o/cm 

1.56 X 10^2 

1.52 X 10®® 

2.37 X 10®° 

to/s 

1.56 X 10® 

6.54 X 10® 

1.42 X 10® 

Tc = c/Cg,0 

5.99 X 10® 

6.62 X 10® 

6.54 X 10® 

TO 

9.12 X 10"^ 

6.41 X 10® 

166.5 

Pr,o/Po 

3.96 

26.32 

13.22 

L/x,y 1 Hq 

0.89 

0.46 

1.17 

Lx/Ho 

6.26 

3.71 

4.60 

Nx,y 

128 

128 

128 

Nx 

768 

1024 

512 


Note: The gravitational acceleration g and radiation flux Fr,i 
coming from the bottom of the simulation box are fixed. The 
optical depth tq, temperature Tq, density poi radiation and gas 
pressure Pr,o-, Pq ^re the fiducial parameters for each run. They 
are also the initial conditions around the iron opacity peak. The 
number of cells are Nx , Ny and Nz . 

copied from the last active zones to the ghost zones used 
for boundary conditions. The vertical component of ve¬ 
locity in the ghost zones is set to be the opposite of the 
velocity in the last active zones. The radiation flux is 
fixed at the constant value F^-^i while the radiation en¬ 
ergy density is integrated to the ghost zones using the 
diffusion equation. For the top boundary, all gas quan¬ 
tities and the radiation flux Fr are copied from the last 
active zones to the ghost zones. When the photosphere 
is included in the simulation box (StarTop & StarMid), 
the radiation energy density in the ghost zones is also 
copied from the last active zones. When the whole sim¬ 
ulation box is optically thick (StarDeep), Er is extended 
to the ghost zones using the diffusion equation. Peri¬ 
odic boundary conditions are used for all quantities in 
the horizontal direction. 

For the short characteristic module we use to calculate 
the VET, the incoming specific intensities from the bot¬ 
tom of the domain are chosen such that angular quadra¬ 
tures of the intensities give a constant vertical flux Fr^i 
and the same Er as in the moment equations at the bot¬ 
tom. The incoming specific intensities at the top of the 
domain are set to zero when the photosphere is inside the 
simulation box. By contrast, when the whole simulation 
box is optically thick, the incoming specific intensities 
from the top of the domain are set to be the same as the 
outgoing specific intensities. 

4. RESULTS 

4.1. The run StarDeep when tq ^ 

The run StarDeep is in the regime where the gas and 
radiation are tightly coupled and we expect efficient con¬ 
vection. Below, we first show how the simulated envelope 
evolves, then summarize its time averaged structure, and 
finally compare the energy transport in our 3D radiation 
hydrodynamic simulation with mixing length theory. 

4.1.1. Simulation History 

Three snapshots of density in Figure]^ show the evolu¬ 
tion of the envelope. The initial density varies only verti¬ 
cally. At the bottom where the temperature is 30% larger 



Figure 3. Snapshots of density for the run StarDeep at time 0 
(left), 21.4fo (middle) and 97.8to (right). The box sizes labeled in 
the plots are in units of solar radius Rq . The horizontal box size 
is Lx = Ly = 20 Rq. Density is in units of po, which is given in 
Table [T] for this run. 


than Tq, the radiation acceleration is smaller than the 
gravitational acceleration and both density and temper¬ 
ature decrease with height. Around z = 13Oi?0, where 
we enter the temperature range of the iron opacity peak, 
the radiation acceleration becomes larger than the grav¬ 
itational acceleration for the fixed radiation flux 
and the density increases with height, peaking around 
z = 16Oi?0. Because the temperature always decreases 
with height according to the diffusion equation, the ra¬ 
diation acceleration eventually decreases after the iron 
opacity peak. When the radiation flux again becomes 
sub-Eddington, the density drops with height. The end 
result is that the initial condition has the density inver¬ 
sion associated with hydr ostatic envelopes with a super¬ 


Eddington radiation flux ( 

Joss et al. 

1973bI. This profile 

is our initial condition for the simul 

ation, and does not 


represent the profile in the MESA ID calculation. 

The middle panel of Figure shows that the initial 
density inversion is unstable. By t « 20to, there are 
large rising and sinking plumes characteristic of convec¬ 
tion and/or the Rayleigh-Taylor instability. We now 
show how this instability can be quantitatively under¬ 
stood as arising from convection in the radiation domi¬ 
nated regime. 

The gas and radiation entropy per unit mass 5'g, Sr are 
given by 


r, _ ks 


In 


Sr = 


AEr 

3^’ 


r p/Po 1 

.ip/po)\ 


and 


(9) 


In our simulations, the total entropy is dominated by the 
radiation. This increases with height near the bottom of 
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Figure 4. Time evolution of the density power spectrum kf(p) 
for the run StarDeep at height 2 = I6O-R0 during the linear growth 
phase of convection. The black and red lines are for the modes 
with wavelength 0.41i^o and O.lSi^o? respectively. The dashed 
black line is a fit of an exponential function to the solid black line, 
which indicates an e-folding time 3.85to for this mode. The power 
spectrum is normalized by the mean density at 2 = 16Oi?0, such 
that the integral of over all wavenumbers is unity. 



Figure 5. History of the radiation energy density weighted 
{Vz^E-r’’ solid line) and density weighted dashed line) ver¬ 

tical velocity for the run StarDeep. 


the simulation box in the initial condition, but necessar¬ 
ily decreases with height near the density inversion where 
T decreases but p increases (see Fig. discussed below). 
We thus expect that our initial condition is convectively 


stable near the bottom of the domain, but convectively 
unstable near the density inversion region in 3D. This is 
exactly what we find in the simulation (see Fig. 1^. Con¬ 
vection happens around 160i?Q within « 20to (equation 
and mixes the initial density inversion. High den¬ 
sity fingers sink and low density gas rises upwards. By 
t = 97.8to (third panel of Figure!^, the bulk of the stel¬ 
lar envelope has reached a statistically steady turbulent 
state driven by convection. 

The initial r adiation Brunt-Vaisala frequency \Nr\ (e.g. 
equation 52 of|Blaes & Socrates|2003|) at the iron opacity 
peak is 0.8/to, but h takes « 2Uio tor convection to begin 
destroying the initial density inversion. The wavelength 
of the fastest growing initial mode is about half of the 
horizontal box size, which also differs from the normal 


Rayleigh-Ta 

ylor instability with a density discontinuity 

(Jiang et al. 

2013a 

) where the shortest wayelengths grow 


stability, we first calculate the 2D discrete Fourier trans¬ 
form in the x, y plane of the density at height z = IGOi?© 
as a function of time. The 2D Fourier coefficient is 
then binned into a ID power spectrum kf{p), where k 
is the magnitude of the horizontal wavenumber. The 
time evolution of two modes with horizontal wavelengths 
O.dliJo = 9.2i?Q and 0.18i7o = 4.Oi?0 are shown in Fig¬ 
ure 1^ The decrease of kf{p) within to probably origi¬ 
nates from an initial vertical oscillation that is present 
due to the imperfect balance between gravity and pres¬ 
sure gradients in the initial condition. After this point, 
the long wavelength mode grows exponentially with an 
e-folding time 3.85to until 15fo- The shorter wavelength 
mode also grows exponentially, but at a rate that is 
smaller by a factor of 3.44 initially. However, it then 
grows quickly after 15fo as convection becomes nonlin¬ 
ear and power from long wavelength modes cascades to 
the short wavelength modes. 

The measured growth rates agree very well with the 
short wa.velengt h WKB linear analysis done by |Blaes fc 


Socrates (2003|) (the last factor of their equation 59). 


PluggiM in numbers for the long wavelength mode in 
Figure and assuming that the vertical wavelength of 
this m^e is about equal to its horizontal wavelength, 
we get a growth rate of 0.27/to, very close to the value 
measured in the simulation. This decline should only 
be valid for wavelengths with growth rates less than the 
magnitude of the Brunt-Vaisala frequency \Nr\, giving 
an expected trans ition wavelength as defined in Blaes| 


& Socrates (2003) of roughly 0.5i7o for this simula¬ 
tion. Hence we expect the longest wavelength convective 
modes to be only weakly affected by radiative damp¬ 
ing, which explains why the measured long wavelength 
growth rate of 0.26/to is of the same order as \Nr\. 
Shorter wavelength modes will be more strongly affected, 
and grow much more slowly. The measured growth rate 
of the shorter wavelength mode here is consistent with 
a reduction factor oc This is consistent with the 

expected behavior given that we are not too far 

from the transition wavelength, and that the short wave¬ 
length mode could have a more horizontal wave vector 
orientation. 

The envelope of the StarDeep simulation shows signif¬ 
icant vertical oscillatory motion. We calculate two mea¬ 
sures of the time-dependent vertical velocity, weighted 
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by energy and by density, 


T/ _ / '^zErdV _ / v^pdV 
“ / ErdV pdV 


( 10 ) 


where the integrals are done over the whole simulation 
box. Histories of these two vertical velocity measures are 
shown in Figure!^ The averaged density, or equivalently 
the total mass, drops by 9% around 20to- At the same 
time, 14,B,, and 14,p reach their peak values. This corre¬ 
sponds to the time when the initial density inversion is 
destroyed due to convection. After the initial transient, 
the envelope slowly settles down to a steady state, and 
both 14 ,Br s-nd Hz,p show oscillations with amplitudes of 
1 — 2% of the radiation sound speed. The average den¬ 
sity also slowly decreases with time because of the mass 
loss through our open top boundary. We have run the 
simulation for more than five thermal times to reach the 
thermal equilibrium. 


4.1.2. Time Averaged Vertical Structures 



100 150 200 


z /^0 


Figure 6. Top: horizontally and time averaged advection flux 
Fadv = '^zEr (black line) for the run StarDeep. T he r ed line is 
the convection flux calculated according to equation (|17| t based on 
the horizontally and time averaged vertical profiles ot temperature 
and pressure in the same run with mixing length parameter a = 
0.55. Mid dle: horizontally and time averaged turbulent velocity V 
(equation [T^, scaled with radiation sound speed Cs (red line) and 
isothermal sound speed Cg (black line) at each height. Bottom: 
vertical profile of V = (V — Vad)/Vad foi" tbe run StarDeep. The 
actual gradient V = d In T/d ln(P + Pr) is calculated based on the 
horizontally and time averaged temperature and pressure profiles. 
The adiabatic gradient is close to 0.25 for the radiation dominated 
envelope. The fact that V 1 implies that the convection is 
almost adiabatic. 


In the 3D radiation hydrodynamic simulations, energy 
can be transported vertically by either photon diffusion 
or turbulent transport. We quantify the turbulent trans¬ 
port using the advection ffux -Fadv, which is 

( 11 ) 


In principle, photon advection can be due to either a 
mean flow (e.g., an outflow) or turbulence. In our simu¬ 
lations, the latter dominates. 

The solid line in the top panel of Figure shows the 
time averaged advective radiation flux J 4 dv a function 
of height where the time average is done between 31to 
and 107to- Most of the flux at ^ < lOOi?© is carried 
by radiative diffusion. In the turbulent region between 
z = llOi ?0 and 1607?©, however, the averaged advection 
flux is more than 40% of Fr,i. We compare this to MLT 
in the next section. 

The turbulent velocity in the envelope can be quanti¬ 
fied as 


V = 


IJ pv■vdV 
JpdV ■ 


( 12 ) 


The middle panel of Figure shows the time averaged 
vertical profile of this turbulent velocity relative to both 
the isothermal sound speed Cg = \/~Pjp and the radiation 
sound speed Cg = ■ The averaged turbulent 

velocity is only 6 — 20 % of the radiation sound speed in 
the turbulent region. Even compared with the isother¬ 
mal sound speed alone, Vjcg is smaller than 0.4. This 
implies that in the efficient convection regime, the tur¬ 
bulent kinetic energy density is much smaller than the 
thermal energy density (< 2 %), as expected. 

Figure shows the quasi-steady envelope structure 
compared with the hydrostatic initial condition. Note 
that the density and temperature are relatively unaf¬ 
fected near the bottom of the box indicating that we 
have placed the bottom boundary deep enough. In the 
final quasi-steady envelope, the iron opacity peak moves 
slightly deeper relative to the initial condition, but the 
density inversion present in the initial condition is gone. 
This is consistent with the ID MESA models which also 
find that for the conditions in StarDeep mixing length 
theory predicts that convection is efficient and is able 
to remove the density inversion. Note that the location 
where the advection flux is largest (z « 12 O 7?0 in Fig. 
™ corresponds to the location of the iron opacity peak in 
Hie final state. This is also where the radiation entropy 
Sr decreases with height, consistent with convection con¬ 
tinuing to drive the turbulent motions and energy trans¬ 
port. As expected, the radiation entropy profile is also 
flatter in the turbulent state after the density inversion 
is removed. 

Because part of the energy flux in the turbulent state 
is carried by advection, the averaged diffusive flux also 
drops below the initial value 74,^. The significant contri¬ 
bution of the advection flux reduces the radiation accel¬ 
eration so that it is smaller than the gravitational accel¬ 
eration. The accelerations due to gas Og and radiation 
Or can be calculated as 


1 dP 

ar = —74,02- (13) 

c 

Notice that only the diffusive radiation flux 74,oz con¬ 
tributes to the radiation acceleration, which is the radia¬ 
tion pressure gradient in the optically thick regime. The 
last panel of Figure [7] shows that gravity is roughly bal- 


74dv — VzEr- 
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Figure 7. Horizontally and time averaged vertical profiles of various quantities for the run StarDeep. Panel (a): density (solid line at top), 
gas temperature (solid black line at bottom) and radiation temperature (solid red line at bottom). The gas and radiation temperature 
profiles overlap as they are strongly thermally coupled. Panel (b): gas entropy Sg (solid black line at top) and radiation entropy Sr (solid 
black line at bottom). These entropies are calculated according to equation §, and both are in units of ke/Ai. Panel (c): total opacity 
ill unit of l/(pof^o) (the solid line). The dashed black lines in these three panels are the initial conditions. Panel (d): volume averaged 
(solid black line) and density weighted (dashed black line) radiation accelerations (ur, flr) as well as the acceleration due to the gas pressure 
gradient {ag, green line). The red line is the sum of dr and ag. In this run with ro ^ Tc, the density inversion in the initial condition is 
gone. 
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anced by the sum of radiation and gas pressure gradients, 
with Gr being at most 90% of g. The fact that < g is 
consistent with the absence of a density inversion in the 
turbulent state. 


4.1.3. Time Averaged Horizontal Structures 

Figure shows horizontal slices of density p and the 
vertical component of Eulerian radiation flux at 
z = 14Oi?0 and time 97.8to for the run StarDeep. There 
is a clear anti-correlation between fluctuations of p and 
Fr,z- When p is larger (smaller) than the horizontally 
averaged density, Fr^z is negative (positive). This is a 
natural outcome of convection, where the Eulerian radi¬ 
ation flux Fr^z is dominated by the advection part VzE^ 
at each location in the regime tq S> Tc- Figure ^ thus 
shows that low density regions buoyantly rise while high 
density fluid elements fall. 

To quantify the level of density and temperature fluc¬ 
tuations, the top panel of Figure shows the vertical 
profiles of standard deviations of p and T scaled with 
the horizontally averaged mean values at each height at 
time 97.8to- The scaled standard deviations peak at 
15Oi?0, where the iron opacity peak is located at this 
time. Beyond the convective region, fluctuations drop 
with height. The scaled density standard deviation is al¬ 
ways smaller than 10% in this case. Temperature fluctua¬ 
tions (< 0.4%) are much smaller. Because radiation pres¬ 
sure Pr ocT"^, the pressure fluctuations are only < 1.6%. 

The cross correlation between two quantities a and b 
can be quantified by the correlation coefficient Ca,b, nor¬ 
malized by their standard deviations <Ja,cri,. The bottom 
panel of Figure shows the correlation coefficients be¬ 
tween p, T and p, F^^z ■ Consistent with the slice shown in 
Figure]^ there is an anti-correlation between p and Fr^z 
in the convective region due to buoyancy. Temperature 
fluctuations also anti-correlate with density fluctuations, 
which means an anti-correlation between density p and 
radiation pressure Pr- In this very optically thick regime 
To ^ Tc, density fluctuations do not reduce the effective 
radiation acceleration, which can be confir me d by com¬ 
paring the volume averaged (equation 131, and the 
density weighted radiation acceleration 


Gr = 


{prttFr,0z} 

"Ap) 


(14) 


Here the average (•) is done horizontally at each height. 
Panel (d) of Figure shows vertical profiles of the time 
averaged radiation accelerations calculated in these two 
different ways. They are almost identical. This explictly 
demonstrates that there is no reduction of the radiation 
acceleration due to density fluctuations in this regime. 

The typical size of the turbulent structures can be 
estimated by calculating the normalized cross correla¬ 
tion Cp^p between p{x,y) and p{x -\- lp,y -\- Ip), where 
p{x -\-lp,y Ip) is the density shifted horizontally with 
distance Ip. The correlation length Ip can be defined to 
to be the shifted distance when Cp^p drops to zero. For 
the run StarDeep, Ip is typically ^ 0.3 — 0.4iJo, with cor¬ 
responding optical depth at the iron opacity peak com¬ 
parable to cjcsp given in Table 


4.1.4. Comparison with Mixing Length Theory 


In the radiation dominated regime, the radiation ad¬ 
vection flux Eadv from the simulations can be directly 
compared with the convection flux from MET. We define 
the ratio of gas pressure to total pressure (/?), the first 
(Fi) and second (r 2 ) adiabatic indices, and specific heat 
at co nstant pressure (cn) fo r a mixture of gas and radia¬ 
tion (Chandrasekhar 19671, and the logarithmic deriva¬ 
tive of density with respect to temperature at constant 
total pressure to be 

P 


P + Pr 


ri = /3 

r 2 = 


(4-3/3)^(7-1) 

+ /3 + 12(y-l)(l-/3)^ 

(4-3/3)ri 


3(l-/3)ri+/3 

1 Fi 


7 - 1/32 
din p 


Q = - 


dlnT 


[/3 + 12(7-1)(1-/3)] 
_ 4-3^ 

p+p,- 


P 


(15) 


where Pr is the zz component of the radiation pressure 
tensor P,.. The adiabatic temperature gradient with re¬ 
spect to total pressure can be calculated as 


Vad = 


dlnT 


dln(P -I- Pr) 


Fa - 1 


(16) 


ad 


The actual gradient V = dliiT/d{hi{P + Pr)) can be 
calculated based on the mean density and temperature 
profiles shown in Figure!^ If we assume that the mix¬ 
ing length is related to the pressure scale height via the 
mixing parameter a, then the convective flux based on 
non-adiabatic mixing l ength theory includ ing the effects 
of radiative cooling is (Henyey et al.||19^ ) 


Qi/"cp kepT fP + P, 
4^2 fi 


1/2 


0,2 (V - V') 


3/2 


(17) 


assuming that the radiation entropy Sr decreases with 
height. Here the gradien t V' is calculat e d base d on equa¬ 
tions (32) and (42)^ in Henyey et al. (1965) given the 
mean vertical profiles of the envelope. In the adiabatic 
limit, V' is reduced to Vad- Then equation[T7]is the same 
as t he formula given by the adiabatic mixing length the¬ 
ory ( |^[I96^ . 

The convective flux calculated in this way with a = 
0.55 is shown as the red line in the top panel of Figure!^ 
which agrees with the time averaged advection flux in the 
simulation. Note that the value of a required to explain 
the flux in the simulations is also quite close to the ratio 
between the density co rrelat ion length and pressure scale 
height i7o (see Section 4.1.3), which is the definition of a 
in mixing length theory. Filially, Figure shows that the 
difference between the actual gradient and the adiabatic 
gradient V = (V — Vad)/Vad is quite small, consistent 
with the expectations of efficient convection. 


^ Notice that there is a typo in equation (42) of [Henyey et al!| 
l|l965||. The a should be a^. We thank Jeremy Goodman for 
pointing this out. 
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Figure 8. Horizontal slices of density p (left) and vertical component of the radiation flux Fr,z (right) at z = 14O-R0 for the run StarDeep 
at time 97.8to- The box size is in units of Rq while units for p and Fr,z are po and carT^ respectively. Notice the strong anti-correlation 
between p and caused by buoyancy. 



100 150 200 


the statistical properties of the turbulent energy trans¬ 
port in our radiation hydrodynamic simulation with tq ^ 
Tc are reasonably consistent with the expectations of mix¬ 
ing length theory, although the value of a is somewhat 
smaller than traditionally assumed in models of massive 


stars (typical values a = 1.0...1.5, see e.g. Yusof et al. 


2013 


Kohler et al. 2015). This is the first numerical 


conhrmation of the adequacy of MLT in a radiation dom¬ 
inated regime. 


4.2. The run StarMid when tq > Tc 

In this regime, rapid radiative diffusion starts to have 
a non-negligible impact on the properties of convection 
in the radiation dominated stellar envelope (the typical 
thermal time scale near the iron opacity peak is ^ 5 
dynamical times). Compared with the previous run, we 
shall see the effects of the reduced optical depth. Note 
that unlike StarDeep, the photosphere is captured within 
the domain of this simulation. 


4.2.1. Simulation History 


z/Re 


Figure 9. Top: vertical profiles of standard deviations of density 
p and temperature T at time 97.8to for fh® run StarDeep. The 
standard deviations are scaled with the mean density and temper¬ 
ature at each height. Bottom: Vertical profiles of cross correlation 
coefficients between p, T and p, Fr,z at the same time. The anti¬ 
correlation between p and Fr,z in this run is due to the buoyancy 
instead of the “porosity” effect, as Fr,z is dominated by the advec- 
tion flux VzEr. 


Taken together, these comparisons demonstrate that 


The initial condition for StarMid has the iron opacity 
peak and density inversion at 78Rq. The density inver¬ 
sion region is, however, convectively unstable and starts 
to mix at time 9.4to as shown in the left panel of Figure 


10 At time 18.8to, the whole envelope is turbulent (right 


panel of Fig. 10). This temporal evolution is very similar 
to the run StarDeep, but the properties of the turbulence 
are different as we now describe. 

Figure shows the radiation energy density and den¬ 
sity weighted vertical velocities {Vz,EzyVz,p) as a func¬ 
tion of time. After the initial density inversion is broken 
around lOfo, 2.4% of the mass is lost through the open 
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Figure 10. Snapshots of density for the run StarMid at time 9.4to 
(left) and 18.8to (right). The box size is in units of Rq while the 
density unit pQ is given in Table ^ The horizontal box sizes are 
Lx = Ly = lOi?0. 



Figure 11. History of the volume averaged vertical velocities 
^z,Er ^rid Vz,p for the run StarMid. 


top boundary. The averaged vertical velocity also reaches 
the peak at the same time. The whole envelope oscillates 
with a period of AIq. The vertical velocities can reach 
10% of the typical radiation sound speed, which is much 


largCT than the oscillation amplitude in StarDeep (Fig¬ 
ure]^. When the envelope expands, it is accelerated by 
the radiation without the density inversion. But because 
of the sensitive dependence of the iron opacity on tem¬ 
perature and the temperature decrease with height, the 
radiation acceleration becomes smaller than g at some 
height and the acceleration stops. The whole envelope 
starts to fall back when the vertical velocity reverses. 



40 60 80 100 
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Figure 12. Top: time and horizontally averaged vertical pro¬ 
files of advection flux (black line) and convection flux predicted by 
mixing length theory (red line). Middle: the averaged turbulent 
velocity V scaled with the radiation sound speed Cs (red line) and 
isothermal sound speed (black line). Bottom: vertical profile of 
the superadiabatic parameter V (black line) and the ratio between 
the local optical depth per pressure scale height t and the critical 
value Tc (red line). These profiles are for the run StarMid and the 
time average is done between 17.6io and 29.2fo. 


4.2.2. Time Averaged Vertical Structure and Turbulent 
Properties 


Figure 12 (solid line) shows that there is significant 


turbulent transport of energy below 65i?0, but the tur¬ 
bulent transport falls off rapidly at larger heights. At 
about the same location 65Rq, the averaged turbulent 
velocity V reaches the isothermal sound speed (middle 
panel of Fig. 12), but is still smaller than the radiation 
sound speed. 

Figure shows the time and horizontally averaged 
vertical pities of density, temperature, opacity, entropy 
and accelerations for this run. The iron opacity peak 
moves inward somewhat in radius, but the density in¬ 
version in the initial condition is largely gone, replaced 
with an extended region below 65i?0 over which the den¬ 
sity changes very slowly with height. Gas and radiation 
are thermally coupled in this region. However, above 
65i?0, the density drops quickly with height and the gas 
temperature starts to deviate from the radiation tem¬ 
perature. Here the local optical depth per pressure scale 
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Figure 13. The time and horizontally averaged vertical profiles of density (top panel of a), temperature (bottom panel of a, solid black 
line for gas temperature and red line for radiation temperature), entropy (panel b), opacity (panel c) as well as accelerations (panel d) for 
the run StarMid. As in Figure!^ the dashed lines in panels (a), (b), (c) are initial conditions while other lines are time averaged profiles. 
The initial condition above 85-R0 is the region above the photosphere where we set it to be isothermal. 


height r becomes smaller than the critical value and 
the photon diffusion time becomes smaller than the lo¬ 
cal dynamic time. This transition to r < Tc is wh y th e 
advection flux drops signifi cant ly in this region (Fig. 12). 

The top panel of Figure [l^ shows the convective max 
calculated based on the time averaged vertical profiles 
and mixing length theory (eq. dzl) with a = 0.45, which 
is again close to the ratio between the correlation length 
Ip and scale height iJo- Below 60Rq where t > Tc, the 
time averaged advection flux is comparable to the mixing 
length theory convection flux, which is consistent with 
Figurej^for the run StarDeep. However, above 65Rq the 


advection flux drops significantly in the simulations, but 
equation ([T7| ) still predicts a very large convection flux 
because or the decreasing radiation entropy (see panel b 
of Fig[T^. This demonstrates that for r < Tc, convection 
becomes inefficient and equation (17) significantly over¬ 
estimates the amount of advection flux with a constant 


a. 


The vertical profile of t/tc, where r = pKtHo is the lo¬ 
cal optical depth per pr essu re scale height, shown in the 
bottom panel of Figure [T^ confirms that the transition 
to inefficient convection happens when t/tc < 1. It also 
happens roughly when the turbulent velocity V becomes 
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Figure 14. Top: vertical profiles of density and temperature stan¬ 
dard deviations scaled with the horizontally averaged mean values. 
Bottom: vertical profiles of cross correlation coefficients between 
p, T and p, FV,z. These profiles are for the run StarMid at time 

18.8to* 


larger tha n th e isothermal sound speed (the middle panel 
of Figure 12), and the superadiabatic parameter V in¬ 
creases sigmncantly to ~ 0.4 (bottom panel of Figure 


12). Notice that throughout the whole envelope, V is 
comparable to Vjcg as in Figure [^for the case StarDeep. 

Panel (d) of Figure shows that the volume aver¬ 
aged radiation acceleration is larger than the gravi¬ 
tational acceleration between 50i?Q and 8 Oi? 0 , which is 
also the region where radiation entropy decreases with 
height. It is thus perhaps surprising that when convec¬ 
tion becomes inefficient above 65i?0 and there is little 
turbulent transport of energy, the envelope does not de¬ 
velop a density inversion again. This is becau se t he den¬ 
sity weighted radiation acceleration dr (eq . 14 sho wn 
as the dashed black line in panel (d) of Figure [T^ is 
significantly smaller than the volume averaged radiation 
acceleration above 65i?0. The sum dr + Ug is compa¬ 
rable to g at the iron opacity peak around 6 Oi? 0 . Beyond 
65i?0, it drops below g significantly. 

Figure clarifies the origin of the reduced den¬ 
sity weighted radiation acceleration. In particular, it 
shows that density fluctuations reduce the effective radi¬ 
ation acceleration above 65i?0 because there is an anti¬ 
correlation between fluctuations of density and radia¬ 
tion flux. The correlation coefficient Cp^Pr ^ in Figure 
[T4| shows this explicitly and is always negative. How¬ 
ever, the nature of the anti-correlations changes at dif¬ 
ferent heights. Below 65i?0, when r > Tc, the situation 
is similar to the run StarDeep, in that buoyancy causes 
the anti-correlation between p and as F!r,z is domi¬ 
nated by the advection flux. Above 65i?0, however, when 
T < Tc, the situation changes. Convection becomes in¬ 
efficient and Fr^z is dominated by the dif fusi ve flux (the 
advection flux drops as shown in Figure 12). The anti¬ 


correlation between p and Fr^z above 65i?0 is thus not 
caused by buoyancy. Instead, it is because the “porosity” 
of the envelope created by the traveling strong shocks al- 
lows radiation to propagate through low density regions 
(|Shaviv|1998|). This reduces the effective radiation accel- 
eration significantly (see panel (d) of Figure 13). Note 
also that in StarMid, the standard deviations of density 
and temperature (scaled with the horizontally averaged 
quantities) are significantly larger than the correspond¬ 
ing values in the run StarDeep (Figure]^. 


4.3. The run StarTop when tq <j; Tc 

This is the regime when gas and photons are loosely 
coupled through the whole iron opacity peak region. 
Convection will be inefficient as in the top region of 
StarMid. The photosphere is now closer to the region 
with the iron opacity peak and accurate radiation trans¬ 
fer is critical. For this run, we use a density floor of 
10“^°/9o to avoid numerical difficulties. 



Figure 15. Time evolution of the density power spectrum kf{p) 
for the run StarTop at height 2: = 13.7i?Q during the first 50fo. 
The solid black and red lines are for the modes with wavelength 
Hq and 0.26i7o respectively. The long wavelength mode has an 
e-folding time 6.25to as indicated by the dashed black line. 


4.3.1. Simulation History 

The iron opacity peak and density inversion are located 
at I3.7i?0 initially. This is again convectively unstable 
and the initial structure is completely destroyed by 50to- 
Just as we found in the simulation StarDeep, this time 
scale is significantly longer than the buoyancy time scale 
I/|A^r|, which is « at the iron opacity peak. We there¬ 
fore again explored the linear growth phase for this run 
by calculating the density power spectrum at the loca¬ 
tion of the initial density inversion. The time evolutions 
of the power kf{p) in two modes with horizo nta l wave¬ 
lengths Hq and 0.26iFo are shown in Figure 15 As in 
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Figure 16. Snapshots of density for the run StarTop at time 49. Ito 
(left) and 138.9to (right). 



Figure 17. History of the volume averaged vertical velocity 14, 
(solid line at bottom) and Vz,p (dashed line at bottom) for the run 
StarTop. The volume average is done for the whole simulation box. 


StarDeep (Fig. |^, the long wavelength mode grows ex¬ 
ponentially from close to the beginning, with an e-folding 
time of 6.25to. However, in contrast to the StarDeep case, 
the short wavelength mode is damped initially and only 
starts to grow quickly after 20to when the e-folding time 
is 3.33to- This later rapid growth is almost certainly due 
to a nonlinear cascade. 

We have examined the behavior of other horizontal 
wavelengths as well, and found that all the modes with 
horizontal wavelengths larger than « OAHq have similar 
growth rates after 12to, while the modes with horizontal 
wavelengths smaller than « 0.3Hq are damped until the 
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Figure 18. Top: time and horizontally averaged vertical profiles 
of advection flux (black line), and convecti on fl ux accordin g to the 
non-adiabatic MLT with a = 0.4 (equation |17|l as given by|Henyey| 


et al. 119651 (red line). Middle: the averaged turbulent velocity V 
scaled with the radiation sound speed Cs (red line) and isothermal 
sound speed (black line). Bottom: vertical profile of superadiabatic 
parameter V. The time average is done between 55.6to and 171.Oto- 
These profiles are for the run StarTop. 


nonlinear cascade sets in. Based on the analysis of Blaes 


& Socrates (2003), we estimate that modes with wave- 
lengths less than about 2.5iJo are in the regime where 
radiative diffusion reduces the growth rate of the convec¬ 
tive instability. This is larger than the size of the box, so 
that all convective modes in the simulation are heavily 



tal wavelength (which is also comparable to the vertical 
width of the density inversion). This agrees very well 
with the measured growth rate. However, shorter wave¬ 
length modes should have smaller growth rates oc 
in contrast to the behavior measured in the simulations. 
It is not clear what is causing this discrepancy, although 
we suspect that the damping of modes with wavelength 
0.3iJo is due to the numerical damping rate being 
larger than the very small growth rate for these modes. 
Interestingly, traveling acoustic waves should also be un¬ 


stable according to equation (62) of Blaes & Socrates 
(2003), with comparable growth rates to the buoyancy 
instability (such waves should not be unstable in the 
StarDeep simulation). Vertically propagating acoustic 
waves (with infinite horizontal wavelength) should grow 
fastest, but these would not necessarily show up as grow¬ 
ing perturbations at a fixed height. We have not further 
investigated this possibility, as it is clear that the unsta¬ 
ble buoyancy modes dominate the evolution observed in 
the simu latio n. 

Figure [T6| shows snapshots of the density structure of 
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the envelope at time 49.1to and 138.9to while the his¬ 
tories of t he vertical velocities are shown in 

Figure The whole envelope becomes turbulent af¬ 
ter the initial density inversion is destroyed around 50to ■ 
Unlike the simulations StarDeep and StarMid, the tur¬ 
bulent envelope is not a relatively quiescent convective 
structure, but instead shows violent, irregular, large am¬ 
plitude oscillations. The oscillation period varies from 
^ 5 to lOto- The amplitude of 14 ,b,. is 0.05cs,o, which is 
already « 18% of the isothermal sound speed. At time 
138.9to when 14,reaches a minimum, the right panel 
of Figure shows that there is even a low density hole 
extending over one pressure scale height Hq. The mass 
loss through the open top boundary is negligible for this 
run. 

According to equation (|^, the typical thermal time 
scale is 4 0.09to for this run. However, for this tq Tc 

regime, 4 is no longer the most relevant time scale for 
some of the dynamics. Instead, the radiation flux be¬ 
haves li ke a roughly consta nt force reducing the effective 
gravity (Jiang et al. 2013a). The more relevant time scale 
is the enective dynamic time scale 


to = 


Og,o 

IffI ’ 


(18) 


where Cg,o is the isothermal sound speed and g = g —dr is 
the effective acceleration due to gravity and the radiation 
force. For an effective acceleration g = —0.06g according 
to the dr shown in panel (d) of Figure [T^ (discussed be¬ 
low), the effective dynamic time scale is tg = 4.6to, which 
is consist ent with the oscillation period of the envelope 
in Figure Note that the effective scale height is then 
Hq = Cg q/I^I = 1.267?0j which happens to be similar to 
Ho- 

4.3.2. Time Averaged Vertical Structures 


Figures 18 and show vertical profiles of various 
time and horizontally averaged quantities in the non¬ 
linear state after 55.64- Unlike the previous two runs 
StarDeep and StarMid, the averaged advection flux is less 
than 0.5% of in the whole envelope as shown in the 
top panel of Figure This is because tq ^ T(. every¬ 
where and gas is not able to trap photons. The averaged 
radiation flux Fr^z is almost a constant vertically ex¬ 
cept very near the top of the domain; the radiation flux 
is nearly constant in time as well, with < 8% variation 
during the whole simulation. 

Because the radiation flux dominates over the advec¬ 
tion flux, the full super-Eddington flux contributes to the 
radiation acceleration of the gas. However, because of 
the large density fluctuations, most photons go through 
regions that have below average density. This red uce s 
the density weighted radiation acceleration dr (eq. 14), 
as shown in panel (d) of Figure 19 H owever, unlik e the 
assumptions made in previous moaels ( Shaviv|1998 ), the 
“porosity” of the envelope due to convection cannot re¬ 
duce the effective radiation acceleration to a value below 
g. At z « 13.3i?0 where the iron opacity peak is located 
now, dr is still about 6% larger than g, while is about 
10% larger than g. This is why the average density profile 
still shows a mild inversion n ear the iron opacity peak, as 
shown in panel (a) of Figure]^ We stress that although 
the density inversion remains, the resulting structure is 


nonetheless very different from the hydrostatic density 
inversion in the initial condition. The high density region 
near the iron opacity peak is continuously reformed and 
then mixed, triggering violent oscillations, shocks, and 
large density fluctuations in the stellar envelope. Physi¬ 
cally, in the regime tq ^ Tc, photon diffusion is so rapid 
that radiation pressure does not respond to the density 
and velocity fluctuations. As the turbulent velocity be¬ 
comes larger than the isothermal sound speed, shocks 
are formed in the envelope (Figure [I^, which causes the 
large density fluctuations. This phenomenon is also ob¬ 
served in radiation magneto- hydrodynamic simu l ations 
of black hol e accretion disks (Turner et al. 2003 Jiang 
et al.|[20l3^ . 


Figure |19| shows that the radiation and gas entropies 


Sr and Sg decrease rapidly with height around the iron 
opacity peak. If we still adopt L/Ho ~ 0.4 for the mixing 
length parameter a, equation ( [l7| ) over-predicts the flux 
relative to the simulations by a factor of more than ^ 3, 
and the location of the peak in the predicted flux also off¬ 
sets from the peak of time averaged advection flux given 
by the simulations, highlighting that existing models of 
inefficient convection in the radiation pressure dominated 
regime are not accurate compared to our simulation re¬ 
sults. The failure of convection to carry a significant 
energy flux is intimately connected with the turbulent 
velocities becoming supersonic, at least relative to the 
gas (isothermal) sound speed. This is expected on theo¬ 
retical grounds and is also consistent with the sim¬ 
ulation results, as shown in the middle panel of Figure 


18 The superadiabatic parameter V is also comparable 
to V/cs and similar to the top region in StarMid (Figure 

T^. 


4.3.3. Properties of the Turbulent Structures 

Figure [^quantifies the large density fluctuations and 
strong anti-correlation between density and radiation 
flux in simulation StarTop. The correlation coefficient 
Cp^Fr^ is almost —1 in the turbulent region between 
12.8i?0 and 13.6i?0, while C'p,T only fluctuates around 
0. This demonstrates that unlike for StarDeep, the anti¬ 
correlation between density and radiation flux is not due 
to buoyancy in the radiation pressure dominated regime 
but is instead due to the photons preferentially diffusing 
through low density channels. The scaled standard de¬ 
viations Sp/p and 5T/T are also much larger than the 
corresponding values in StarDeep and StarMid, which is 
consistent with the images in Figure [l6| 

The typical size of the turbulent structures can also 
be quantified by the correlation length Ip, which varies 
between OAHq and 0.5Hq. This is likely set primarily 
by geometry as in standard mixing length arguments. In 
particular, the optical depth across the correlation length 
Ti is much smaller than Tc, which means that photon 
diffusion is very rapid for all of the turbulent fluctua¬ 
tions and cannot be important for setting a characteristic 
scale. The correlation length Ip is also much smaller than 
the characteristic turnover wavelength of acoustic waves 
driven unstable by the dependence of opacity on density 
(Blaes & Socrates 20031. Once again, it is simply the fact 
that the convective turbulence is supersonic with respect 
to the isothermal sound speed in the gas, and that radi¬ 
ation pressure cannot respond due to the rapid diffusion. 




















16 


Y.-F. Jiang et al. 



12.5 


13 13.5 

z/Re 


14 






12.5 


13 13.5 

z/Re 


14 



z/Re 


Figure 19. The time and horizontally averaged vertical profiles of density (top panel of a), temperature (bottom panel of a, black and 
red lines for gas and radiation temperature), entropy (panel b), opacity (panel c) as well as accelerations (panel d) for the run StarTop. 
The dashed lines in (a), (b) and (c) are initial conditions while other lines are time averaged profiles between 55.6to and 171.Otp. 


that is producing these nonlinear structures. 

5. DISCUSSIONS AND CONCLUSIONS 

We have presented three-dimensional (3D) radiation 
hydrodynamic simulations of the radiation pressure dom¬ 
inated envelopes of massive stars near the iron opacity 
peak. The resulting dynamics depends sensitively on the 
ratio of the photon diffusion time to the dynamical time. 
This is set by the ratio of the local optical depth tq to 
the critical optical depth Tc = cjcg, where c is the speed 
of light and Cg is the isothermal sound speed. Our sim¬ 
ulations cover the regimes rp (StarDeep), tq ~ Tc 

(StarMid), and rp <C Tc (StarTop). 


The iron opacity peak at T « 1.6 x 10^ K leads 
to locally super-Eddington fluxes in the envelopes of 
massive stars. In radiative equilibrium this produces a 
pronoun ced density inversion in one-dimensio nal stellar 


models (Joss et al. 1973b Paxton et al. 


2013). The re¬ 


sulting density inversions are convectively unstable, but 
the properties of the resulting convection zones have un¬ 
til now been very uncertain when rp < Tc and convection 
is expecte d to be rela tively inefficient. 

Figures and compare the time and horizontally 
averaged turbulent properties of our three runs, show- 
ing the gas, radiation, and kinetic energy densities (Fig. 


21) and the ratio of the convective and diffusive energy 
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Figure 20. Top: vertical profiles of density and temperature stan¬ 
dard deviations scaled with the horizontally averaged mean values. 
Bottom: vertical profiles of cross correlation coefficients between 
p, T and p, FV, 2 - These profiles are for the run StarTop at time 

138.9£o- 


Figure 22. Vertical profiles of the ratio between time averaged 
advection (convection) flux Fadv = '^zEr and diffusive flux 
for the three runs. Convection transports significant energy in 
StarDeep and in the deeper layers of StarMid. But in the outer 
layers of StarMid and throughout StarTop, photon diffusion is so 
rapid that it strongly suppresses the turbulent transport of energy. 




(z-Zo)/H^ 

Figure 21. Time and horizontally averaged vertical profiles of 
radiation energy density (Er-, red lines), gas internal energy density 
{Eg, solid black lines) and kinetic energy density green lines). 
From top to bottom, they are for StarDeep, StarMid and StarTop 
respectively. The black dashed vertical lines in each panel indicate 
the time averaged location of the iron opacity peak. The horizontal 
axis is the distance from the bottom boundary zq in units of the 
scale height Hq for each run. 


fluxes, respectively. 

When To Tc as in StarDeep, convection is efficient 
and the averaged radiation advection flux is consistent 
with the convection flux as calculated using the MET 
formalism (equation [l7|. We find that the mixing length 
parameter required to match the advection flux in the 3D 
simulation is a = 0.55, roughly consistent with the typi¬ 
cal size of the turbulent structures Ip/H q « 0.3 — 0.4 (as 
assumed in mixing length theory). The density inversion, 
present in the ID radiative equilibrium initial condition, 
is absent in the turbulent state. This is because a signifi¬ 
cant fracti on o f the luminosity is transported by convec¬ 
tion (Fig. [2^ , which reduces the radiation acceleration 
to be sub-Eddington. For StarDeep, the turbulent ve¬ 
locity is very subsonic with the turbulent kine tic energy 
less than 1% of the thermal energy (Fig. 21 1 . The dif¬ 
ference between the actual thermal gradient V and the 
adiabatic gradient Vad is smaller than 8 % of Vad- All 
of these properties are consistent with expectations for 
efficient convection. 

When To < Tc, the advection flux drops signihcantly 
as photons cannot be trapped by the fluid due to rapid 
diffusion. We see this transition in the surface layers of 
run StarMid. As shown in Figure 12 below 6 Oi? 0 , the 
aver aged radiation advection flux is consistent with MET 
(eq. 17). The turbulent velocity is subsonic and V ^ 
1 as m StarDeep. However above 607?©, equation (17) 
signihcantly overpredicts the average convection flux if 
we use a constant a, a choice motivated by the fact that 
the correlation length of the turbulence {Ip/H q) is almost 
the same through the whole envelope. Above 60i?Q the 
turbulent velocity reaches the isothermal sound speed 
and V increases signihcantly. 
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The run StarTop has tq <C Tc and convection is very 
inefficient. In this case the average radiation advection 
flux is <C 1% of the diffusive flux throughout the envelope 
(Figu re [2^ . If we adop t the non-adiabatic MLT model 
(e.g.jHenyey et al.|iyb5l|Magic et al.|2015|) with a = 0.4, 
which IS the ratio of the turbulent correlation length to 
the pressure scale height, we find that it still significantly 
over-predicts the convection flux (Figure 181. The differ¬ 
ence is probably caused by a combination or the porosity 
effect and vertical oscillation of the envelope as explained 
blow. 

The turbulent velocity in StarTop is larger than the 
isothermal sound speed and strong shocks are formed, 
driving lar ge density fluctuations. The bottom panel of 
Figure shows, however, that the kinetic energy den¬ 
sity only reaches the gas internal energy density, which 
is still much smaller than the radiation energy density. 
The large density fluctuatio ns allow radia tion to escape 
somewhat more efficiently (Shaviv 1998) than through 
a homogenous medium (“porosity”), reducing the effec¬ 
tive radiation acceleration. This effect is not, however, 
large enough to make the radiation acceleration smaller 
than the gravitational acceleration. The density inver¬ 
sion present in the initial conditions is initially destroyed 
by convection and the envelope expands, but the gas 
quickly falls back because of the reduced opacity and ra¬ 
diation acceleration. The density inversion then reforms. 
This cycle repeats, with the envelope undergoing strong 
vertical oscillations with a period of a few hours. The 
density inversion predicted by ID radiative equilibrium 
models still exists in the time averaged vertical profile of 
StarTop, in contrast to StarDeep and StarMid where con¬ 
vection is able to largely eliminate the density inversion. 

5.1. Observational consequences 

A signature of inefficient convection is the appear¬ 
ance of supersonic turbulent velocities with respect to 
the isothermal sound speed (but not the radiation sound 
speed). This is because rapid diffusion leads the photons 
and gas to decouple and the radiation pressure cannot re¬ 
spond to density fluctuations. This in turn makes the gas 
highly compressible, leading to shocks and large density 
fluctuations. 

In the calculations that include the photosphere (Star- 
Top and StarMid) we observe turbulent velocities reach¬ 
ing the isothermal sound speed close to the stellar sur¬ 
face (« 50 km/s). These large velocities can impact the 
spectroscopic measurements of fine profiles, as quantified 
by the macroturbulence (microturbulence) parameters. 
These parameters quantify the velocity fields at the stel¬ 
lar photosphere with correlation length larger (smaller) 
than the line-forming region. The extent of the fine 
forming region is comparable to the local pressure scale- 
height, which is about the correlation length of the turbu¬ 
lent structures in our calculations. Therefore we expect 
the velocity fields observed in the simulations to poten¬ 
tially impact both micro- and macro- turbulence mea¬ 
surements. This is important, as the broadening of spec¬ 
tral fines is used to measure the (projected) rotational ve¬ 
locity of stars, which in the presence of extra turbulent 
broadening could be be systematically higher tha n the 
actual rotation rate ( Ramirez-Agudelo et al.pOlS ). The 
idea that inefficient convection at the iron opacity peak 
could be responsible for the extra broadening of spectral 


lines in hot massive s tars has been discussed previously 
( Cantiello et al.|[^009[ ), where they focused on the role of 
gravity waves excited by turbulent convection. Recently 
the potentially observable signature of turbulent pres¬ 
sure in the inefficient outer con vective regions of massiv e 
stars has been highlighted by Grassitelli et al. (2015), 
who showed a correlation between observed macroturbu- 
lence and turbulent pressure (as calculated in the context 
of ID MLT). 

The vertical oscillations of the envelope in StarMid and 
StarTop also cause the radiation flux coming from the 
photosphere to vary with time. The standard deviations 
of the total radiation flux leaving from the photosphere 
are « 13% and « 2.7% of the time averaged mean values 
in StarMid and StarTop respectively. This means that 
this oscillation can be observed in principle. However, 
direct comparisons with observations require global cal¬ 
culations of the envelope. 

Strange mode instabilities are commonly associated 
with opacity bumps and d ensity inversions in ID mod¬ 


els o f massive stars ( e.g., Kiriakidis et al. 1993 Saio 


2009 Saio et al. 2013). These standing wave instabili¬ 


ties are fundamentally acoustic in nature, and are gener¬ 
ally associated with a trapping region around the opacity 
bump/density inversion. The excitation mechanism of 
these modes is due to opacity variations in the acous¬ 
tic wave, and is related to the growth mechanism of 
traveling acoustic waves discussed by Blaes & Socrates 
(2003). As we briefly noted above in section 4.3.1, con¬ 
ditions in StarTop at least are such that these traveling 
wave instabilities should be driven, and the fact that the 
density inversion persists in the time-averaged structure 
may permit the existence of trapped strange modes. It 
would be interesting to do a global stability analysis on 
the time-averaged structure of StarMid and StarDeep, to 
see if standing acoustic wave instabilities could still be 
trapped even in the absence of density inversions. 

5.2. Implications for ID stellar evolution 

For To > Tc, we find that mixing length theory provides 
a reasonable description of the convective energy trans¬ 
port in the 3D simulations. However, the mixing length 
parameters we find {a « 0.5) are noticeably smalle r than 
the canonically adopted valu es {a > 1, see e.g. Yusof 


et al. 2013 Kohler et al.|2015|). We also find that as tq/Tc 
decreases, the mixing length parameter decreases. This 
is exactly the opposite of what has been proposed by dif¬ 
ferent groups dealing with the numerical difficulties asso¬ 
ciated with radiation dominated stellar envelopes in ID. 
They suggested that convection in this situation should 
be artificially enhanced (equivalent to increasing the a 
parameter with respect to the canonical values above), 
justifying this choice with eithe r the need to remove an 
“unphysica l density inversion” (Maeder 1987| Ekstr6m| 
et al.|2012 ) or the need to incl ude some (unknown ) extra 
energy transport mechanism ( Paxton et al.|[2M3 ). 

Moreover we find that density inversions are present in 
the time averaged structure of the stellar envelope when 
To ^ Tc, i.e. when an opacity peak is close to the stellar 
surface. In the case of the iron opacity peak, this occurs 
for very massive main sequence stars, though we expect 
the same physics to occur in cooler objects due to the 
effect of the H and He opacity bumps (not studied in 
this work). The behavior of these density inversions is 
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quite different as compared to the steady ID solution, as 
they tend to be cyclically destroyed and reformed and 
are associated with very large density inhomogeneities. 
One important effect that is currently missing in the ID 
models is the inclusion of a “porosity factor”, reducing 
the effective rad iation acceleration when tq < (e.g.. 


Paczyhski 1969). The dependency of the porosity on 
optical depth, resolution and magnetic fields will be the 
focus of future work. 

The turbulent, non-stationary behavior of radiation 
dominated envelopes suggests a possible modification of 
the stellar mass loss rel ative to spherically symmetric 
hydrostatic models (e.g., Owocki et al. 2004). Unfortu¬ 
nately as we only simulate a local patch of the envelope, 
we cannot determine the fate of the mass that leaves our 
simulation box. The mass loss rate also cannot be calcu¬ 
lated self-consistently based on these simulations. This 
requires global calculations of the envelope to capture 
the change of gravitational acceleration and geometric 
dilution correctly. 

The majority of mas sive stars are found in binary sys¬ 
tems (Sana et al.||2012). Depending on their orbital sep- 
aration and the evolution of their radii, these stars can 
fill their Roche lobe and strongly interact alr eady during 
the main sequence (e.g., de Mink et al.|[20l4 ). As shown 
in this paper, the energy transport in the radiation dom¬ 
inated envelopes of these stars is not properly modeled in 
current ID calculations, resulting in very uncertain val¬ 
ues for the stellar radii. This could have important im¬ 
plications for the evolution of massive binaries and, for 
example, the predicted rates of BH-BH binary me rges, 


which are promisin g sources of gravitational waves (Bel- 
if t 


czynski et al. 


2014b 


One limitation of the simulations presented here is that 
we have neglected the effects of the stellar magnetic field. 
Recently, magnetic buoyancy has been found to signifi¬ 
cantly enhance the radiation ad vection flux, even with a 
relatively we ak magnetic field ( jBlaes et al. 2011 Jiang 


et al. 2014a). A fracti on of massive star s do show sig- 


nihcant magnetization (|Wade et al. 2014), and dynamo 
generated magn etic helds could be present in the iron 
convection zone (|Cantiello & Braithwaite||2011|). There¬ 
fore the effects of magnetic field should be studied in 
future work. 
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